/*-------------------------------------------------------------------------------*
 *              Copyright (C) 1998 Texas Instruments Incorporated.               *
 *                              All Rights Reserved                              *
 *-------------------------------------------------------------------------------*/
/*===============================================================================*
 * NOTES:  The filter coefficients generated with this algorithm can be used for *
 *         scale_horz() and scale_vert().                                        *
 *===============================================================================*/
#include "filt_gen.h"

#define W  (0x800)

/* 
    p / q          scale factor
    4              number of taps in filter
    guide_table    ptr to the table of offsets along the filter array
    filter         ptr to the resulting array of filters this has dimension 
                   p * L
*/

void filter_design(int p, int q, short * guide_table, double * filter)
{
    int i,j,k;
    double tap, fp, theta;
    j = 0;
    
    fp = (double)p;
    for (i = 0; i < p; i++)
    {   
       for (k = -1; k < 3; k++)
       {
           theta = (double)(q*i - p*k - j*p );
           theta = theta + 0.5*fp;
           theta = theta/fp;
           if (theta < 0.0) theta = -theta;
           tap = 0.0;
           if (theta < 2.0)
           {
               tap = -0.3889*theta*theta*theta + 
                      2.0000*theta*theta -
                      3.3333*theta + 
                      1.7778;
           }
           if (theta < 1.0)
           {
              tap = 1.1667*theta*theta*theta - 
                    2.0000*theta*theta +
                    0.0000*theta +
                    0.8889;
           }
           *filter++ = tap;
       }
       j += ((short)guide_table[i]);
    }
}

void filter_normalize(int p, double * filter)
{
  int i,j;
  double nom;
  for (i=0; i < p; i++)
  {
    nom = 0.0;
    for (j=0; j < 4; j++) nom += filter[i*4 + j];
    if (nom < 0.0) nom = -nom;
    for (j=0; j < 4; j++) filter[i*4 + j] = filter[i*4+j]/nom;
  }
}

void bresenham_guide(int p, int q, short * guide_table)
{
    int j;
    double stride, i, fp, fq, delta, ratio, step, sum;

    fp = (double) p;
    fq = (double) q;

    j = 0;
    delta = fq/fp ;

    /* for 27/100 you want 4 ~ 100/27 */
    /* for 4/3 you want 1 ~ 3/4 */

    step = (int) (1 + delta);
    sum = step;
    for (i = 1.0; i < p+1.0; i = i + 1.0)
    {
      ratio = sum/i;
      if (ratio > delta)
      {
        stride =  step - 1.0;
        guide_table[j] = (short) stride;
        sum += stride;
      } else {
        stride =  step ;
        guide_table[j] = (short) step;
        sum += step;
      }
      j++;
    }
}


void unrollx1 (short * jump0, short * jump1, int t_l_hh, int t_n_hh)
{
    int j, k;
    for (j = 0; j < t_n_hh; j++)
    {
        if (t_l_hh > 4)
            for (k = 4; k < t_l_hh; k += 4) *jump1++ = 4;
        *jump1++ = jump0[j] - t_l_hh + 4;
    }
}

void unrollx2 (short * jump1, short * jump2, int t_l_hh, int t_n_hh)
{
    int i, j, k, l, M;
    short * jump1_e, * jump1_o, * jump2_e, * jump2_o;

    k = (t_l_hh)>>2;
    l = 0;
    M = (t_n_hh * t_l_hh)>>2;
    jump1_e = jump1;
    jump1_o = jump1;
    jump2[0] = 0;
    jump2[1] = 0;
    for (i=0; i < (t_l_hh)>>2; i++) jump2[1] += 2*jump1[i];
    jump2_e = jump2 + 2;
    jump2_o = jump2 + 3;
    for (j=0; j < t_n_hh>>1; j++)
    {
        for (i = 0; i < (t_l_hh - 4)>>2; i++)
        {
            *jump2_e = 2*jump1_e[l]; l = (l+1)%M; jump2_e += 2;
            *jump2_o = 2*jump1_o[k]; k = (k+1)%M; jump2_o += 2;
        }
        *jump2_e = 0;
        *jump2_o = 0;
        for (i=0; i < (t_l_hh + 4)>>2; i++)
        {
            *jump2_e += 2*jump1_e[l]; l = (l+1)%M;
            *jump2_o += 2*jump1_o[k]; k = (k+1)%M;
        }
        jump2_e += 2;
        jump2_o += 2;
    }
    *jump2_e = jump2[2];
    *jump2_o = jump2[3];
}
 
/* the filters are interleaved even, odd, even , odd, etc, in groups of
   4 taps. for each output point in the table. */
/* for l_hh = 4 this is the same as before */
 
/* p/q is the scale factor both must be divisible by 2 */
/* hh is the pointer to the filter bank                */
/* jump_table is the horizonbtal jump table            */
/* guide table is the verticle                         */

void filt_gen(int p, int q, short * hh,
              short * vguide_table, short * hguide_table)
{
 int j, k;
 short round, tap;
 double filter[200*4];
 short  jump_table[200];
 int l_hh = 4;
 int n_hh = p;

 bresenham_guide  (p, q, vguide_table);
 unrollx1 (vguide_table, jump_table, l_hh, n_hh);
 unrollx2 (jump_table, hguide_table, l_hh, n_hh);
 filter_design    (p, q, vguide_table, filter);
 filter_normalize (p, filter);

/* for the horizontal filter to avoid bank conflicts */

 for (k = 0; k < p; k++)
 {
    for (j = 0; j < 4; j++) 
    {
      tap = (short) ((0x4000)*filter[k*4+j]);
      round = 4;
      if (tap < 0) round = -3;
      tap = (tap + round)>>3;
      *hh++ = tap ;
    }
 }

}

/*-----------------------------------------------------------------------------*
 *              Copyright (C) 1998 Texas Instruments Incorporated.             *
 *                              All Rights Reserved                            *
 *----------------------------------------------------------------------------*/
 
